Stability and excitations of a bilayer of strongly correlated dipolar Bosons 
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We study correlation effects and excitations in a dipolar Bose gas bilayer which is modeled by 
a one-dimensional double well trap that determines the width of an individual layer, the distance 
between the two layers, and the height of the barrier between them. For the ground state calculations 
we use the hypernetted-chain Euler Lagrange method and for the calculation of the excitations we 
use the correlated basis function method. We observe instabilities both for wide, well-separated 
layers dominated by intra-layer attraction of the dipoles, and for narrow layers that are close to 
each other dominated by inter-layer attraction. The behavior of the pair distribution function leads 
to the interpretation that the monomer phase becomes unstable when pairing of two dipoles becomes 
energetically favorable between or within layers, respectively. In both cases we observe a tendency 
towards "rotonization" , i.e. the appearance of a soft mode with finite momentum in the excitation 
spectrum. The dynamic structure function is not simply characterized by a single excitation mode, 
but has a non-trivial multi-peak structure that is not captured by the Bijl-Feynman approximation. 
The dipole-dipole interaction between different layers leads to additional damping compared to the 
damping obtained for uncoupled layers. 



I. INTRODUCTION 

Experimental advances in achieving Bose-Einstein con- 
densation (BEC) of atoms with large magnetic moment 
( 52 Cr [JJ |2], 164 Dy g], 168 Er, g]) and in generating 
quantum gases of heteronuclear molecules (KRb [ME], 
LiCs [9] , LiK [10] , RbCs HHH2I) by Feshbach association 
have lead to a growing interest in effects caused by the 
dipole-dipole interaction (DDI) . The shape of a trapped 
dipolar condensate [12] [T3] and its stability against col- 
lapse [2j [15] have been investigated. The dynamics [TBl - 
|2"U] has been studied theoretically, but recently also ex- 
perimentally [2JJ. The generation of novel phases with 
topological order using polar molecules has been pro- 
posed [22 . A recent review of the field can be found 
in Ref. 23, The strength of the DDI can be character- 
ized by the dipolar length r$ — mCdd/i^TtJi ), where m 
is the mass of the dipolar atom or molecule and Cdd is 
proportional to the square of the dipole moment. For 
magnetic moments, r is usually much smaller than for 
electric dipole moments of molecules, which can range 
to thousands of A. Since achieving BEC with heteronu- 
clear molecules is much harder than with homonuclear 
molecules, Er 2 is a promising candiate to reach a much 
stronger DDI regime with purely magnetic dipole mo- 
ments [24] . 

The two-dimensional limit of a dipolar Bose gas (DBG) 
polarized perpendicularly to the plane has been studied 
extensively by quantum Monte Carlo methods [18 ] l25lf28] 
for a wide range of dimensionless densities nr^, includ- 
ing high densities where the dipole-dipole repulsion leads 
to such strong in-plane correlations that the excitation 
spectrum exhibits a roton similar to the roton in super- 
fluid 4 He, and even higher densities where the ground 



state of the 2D DBG is a triangular crystal. Such large 
nr 2 . may soon be in experimental reach because r$ can be 
very large, as mentioned above. For the 2D DBG with 
tilted polarization a stripe phase has been predicted re- 
cently for sufficiently large tilt angle [30J [5S] . Also the 
more complicated case of a quasi-2D layer was studied, 
i.e. of a DBG in a one-dimensional trap U ex t(z). The 
finite extent in this direction allows pairs of particles 
to explore the anisotropy of the DDI. Already by using 
the mean field approach of the Gross-Pitaevskii equa- 
tion, "rotonization" of a quasi-2D layer of a polarized 
DBG was found to occur if the strength of the DDI sur- 
passes a critical value with respect to the short-range 
repulsion [T7] [3D]. The roton in this case is not a signa- 
ture of the repulsive correlations as in the 2D limit for 
high densities, but a signature of the attractive correla- 
tions for head-to-tail configurations of pairs of dipoles. 
We have shown this conclusively in Ref. 1191 where even a 
cross-over between these two kinds of rotons was demon- 
strated by variation of the trap length. While the effect 
of attractive correlations on the excitation spectrum can 
be qualitatively described by mean field methods that 
optimize only the ground state density by the Ritz' varia- 
tional principle, repulsive correlations leading to 4 He-like 
rotons require the optimization of at least density and 
pair density. This is exactly what the family of hyper- 
netted chain Euler-Lagrange (HNC-EL) methods does. 
The HNC-EL method was therefore used in Ref. [TpJ in 
order to investigate both kinds of rotons. A summary 
of how HNC-EL works is given in the next section, all 
details about the method can be found in Ref. [31]. 

In this work we extend our previous investigation of 
a single quasi-2D DBG layer [19l [32] to a bilayer. The 
coupling between layers via the long-ranged DDI and the 
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possibility of pairing of dipoles on different layers to form 
dimers (and more generally of n bound dipoles in n lay- 
ers) has been investigated previously [33J El] • Superflu- 
idity of fermionic bilayers was studied in Ref. [35 . The 
bound and scattering states of just two dipoles on differ- 
ent layers have been studied |36j |37j as well. 

Our bilayer is realized by a one-dimensional double well 
potential 

Uext(Zi) = A {cos (qzi - it) + A cos (2qz t - 2ir)} . (1) 

A DBG in this trap is homogenous and infinitely large in 
x- and ^-direction and finite in the confinement direction 
z. The dipole moments are aligned along the z-direction, 
therefore the dipole-dipole interaction (DDI) potential 
takes the form 



C dd l-3cos 2 ft 
4-7T r„ 3 „. 



(2) 



where ftj is the angle between the dipoles i and j mea- 



sured from the z-axis, and r. 



To stabilize the 



system against collapse [2] we add a hard core repulsion, 
that is modeled by a r"- 1 potential. The Hamiltonian de- 



scribing this many-body system, in the reduced length 
and energy units, vq 
respectively, looks as follows: 



mCddl{^Ti 2 ) and e = /(mrl) 



i=l i=l 



(3) 
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Correspondingly, we give all values for energy, length, 



wave number, and density in units of eo, rp, 



and 



; hence all quantities in tables and figures are dimen- 
sionless. 

In previous work [16j H3 113 E2] it was established 
that a translationally invariant single layer of a DBG 
in a one dimensional harmonic trap can become unsta- 
ble due to the attractive part of the interaction. The 
pair distribution calculated in Ref. [TJJ] shows that this 
instability can be understood as a dimcrization, where 
two dipoles can form a bound state; such weakly bound 
dipoles would not be stable and indeed experiments with 
52 Cr [5] and 168 Er g] observe a collapse of the BEC for 
traps that are too wide in the polarization direction. Sim- 
ilarly, two coupled layers can become unstable not only 
due to the attractive interaction within a layer, but also 
due to the attraction between dipoles in different layers. 
This latter "instability" actually indicates the dimeriza- 
tion of dipoles in different layers. As long as the barrier 
between the layers is high enough that the two bound 
dipoles remain in their respective layer, such a dimerized 
phase would be stable. 

We have studied the DBG in the trap potential ([I ) for 
various potential parameters A, q, and A that contro the 
barrier height between the wells, their separation, and 



their individual width. We changed the parameters such 
that we can study the transition from two broad, but 
well-separated layers to two thin, but close layers. We 
thereby go from a limit that is dominated by mira-layer 
attraction to a limit dominated by mter-layer attraction. 
Both limits are characterized by the appearance of a soft 
mode (roton) with a respective typical parallel wave num- 
ber fcyx = 0(1). In the first limit, x = a,ho is the oscillator 
length of the approximately harmonic well felt by each 
layer; in the second limit, x = d is the distance between 
the layers. The six combinations of potential parameters 
A and q that we used are listed in the first two columns of 
table |llj A was fixed to A = 0.3. The corresponding trap 
potentials are plotted in the lower panel of Fig. [I] where 
we scale it by 100/^4 in order to show all six potentials 
in the same figure. 

We note that, although the average of the DDI over 
the whole plane vanishes, two dipoles on different two- 
dimensional planes separated by a distance d always 
dimerize, i.e. form a weakly bound state due to the at- 
tractive head-to-tail well of the DDI, regardless of the 
value of d [3H]. Hence in the zero density limit, the 
ground state is always dimerized. At finite density, our 
calculations in the 2D limit show that dimerization is 
suppressed by many-body effects [59] . In other words, 
increasing the density in the two layers stabilizes the 
monomer phase. 



II. GROUND STATE 

Ground state properties of Bose gases can be calcu- 
lated using various methods, such as the Gross-Pitaevskii 
method [40-42 and quantum Monte Carlo methods [43T 
|4"5] . The Gross-Pitaevskii method is widely applied for 
dilute systems, where the correlations between particles 
are sufficiently weak such that the interaction between 
them can be approximated by an effective mean potential 
felt by each particle. Quantum Monte Carlo on the other 
hand is also suited for strongly interacting systems, but 
is computationally demanding. For our calculations of 
the ground state properties we use the hypernetted-chain 
Euler Lagrange (HNC-EL) method, which is a variational 
method suitable for strongly correlated systems [3TJ , but 
with lower computational demands than QMC. Starting 
point is a Jastrow-Feenberg ansatz for the many-body 
wavefunction 
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which is optimized by solving the Euler-Lagrange equa- 
tions numerically 
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Here p(r) = pi(r) is the one-body density and g(r±, T2) 



P2(*i, r 2) (p( r i)p( r 2)) 1 is the pair distribution function. 
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0.20 
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4.69 


0.16 
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3.04 





TABLE I. Table containing the potential parameters A and 
q, the corresponding tunnel splitting for a single particle and 
the chemical potential p of the many-body system. 



pi(r) and P2(ri,ra) are special cases of the n-body den- 
sity reduced from the full ./V-body probablity of the (nor- 
malized) full wave function "0o( r ii • ■ ■ j rjv) 

p„(ri, . . . ,r„) = 

(N^ny. J d3r " +1 ■ -- d3r N l^o(ri, . ■ • ,rAr)| 2 

Due to translational invariance in x and y-direction, 
for the present layer geometry all two-body functions, 
such as the pair distribution function, depend on three 
variables: the mo dulus of the projection o f r = ri — r 2 on 
the plane, r\\ = \J [x\ — x 2 ) 2 + (yi — II2) 2 , and the two z- 
componcnts Z\ and z-i- Hence we effectively have a pair 
distribution function g(z\, Z2, rn). 

All calculations are done for a total area density of 
Tiro = 1, i.e. 717-q = 1/2 for each of the two layers. 
The HNC equations can be formulated in terms of Mayer 
cluster diagrams known from classical statistical mechan- 
ics [?7]. The exact solution of the HNC equations would 
require the calculation of a class of diagrams called ele- 
mentary diagrams that cannot be summed exactly. Ele- 
mentary diagrams are especially important at high den- 
sities, while they can be neglected a lower densities. Fur- 
thermore, we note that we use a Jastrow-Feenberg ansatz 
that does not extend to triplet correlations, u^fa, rj, r^). 
From He we know that both elementary diagrams and 
triplet correlations are important for quantitative agree- 
ment with experiment and quantum Monte Carlo simu- 
lations [48]. For the area density used in this work, we 
have checked the influence of the elementary diagrams 
and the triplet correlations in the 2D limit. In this limit 
correlations are stronger than for quasi-2D geometries, 
so the 2D limit gives a conservative estimate of their 
importance. We found that they improve the accuracy 
of the static structure function S(k) (see below) by less 
then 2%. Therefore we neglect elementary diagrams and 
triplet correlations. In Ref. |32j we compared results for a 
2D system of aligned dipoles to QMC calculations [26 for 
the density po = 2. The agreement between the HNC-EL 
results and the QMC calculations for the static structure 
function is very good. The total energy for both calcula- 
tions differs by about 2.8%, if elementary diagrams and 
triplet correlations are neglected. For the even smaller 
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FIG. 1. (Color online) The density profile p(z) is shown in 
the upper panel and the trapping potential U ex t (z) according 
to eq. |T| is shown in the lower panel. The respective trap 
parameters q and A are listed in the upper panel. U e xt(z) was 
scaled by a factor of ^ in order to show all potentials in the 
same figure. 



area density considered here, we expect the accuracy to 
be even better. 

In table [H] we show the six parameter combinations 
that we choose for the trap potential ([T]) along with the 
tunnel splitting for a single particle and the chemical 
potential p of the many-body system, p is measured 
with respect to the single-particle ground state energy, 
i.e. with respect to a non-interacting Bose gas in the 
same trap potential. As expected, the chemical potential 
increases if we decrease the thickness of the individual 
layers, which is due to the intra-layer repulsion of both 
the DDI and the short-range interaction. With increas- 
ing parameter q, the two trap wells are not only getting 
closer, but, with our choice of parameter combinations, 
also the tunnel splitting increases. As mentioned above 
we gradually move from thick layers that are widely sep- 
arated to thin layers that are close to each other, while 
keeping the total area density fixed at nr 2 = 1. The re- 
sults for the density profiles p{z) of the DBG in the trap 
potentials can be seen in the upper panel of Fig. [l] The 
corresponding trap potentials U ex t{z) are shown in the 
lower panel using the same line style and color (online). 
U ex t{z) is scaled by the inverse of the trap parameter A 
such that all six potentials can be shown using the same 
scale. Each potential is offset such that the ground state 
energy of a single particle is zero. 

In the limit of zero density or in the non-interacting 
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FIG. 2. (Color online) The density profile p(z) is shown for 
the three wider traps and compared with the single-particle 
density |</>o(-z)| 2 in the same traps (dotted lines). 



FIG. 3. (Color online) Static structure function S[ku) as a 
function of the parallel momentum k\\ for the six different 
trapping potentials. 



limit, the density profile is given by the square of 
the ground state solution, |0o(z)| 2 , to the one-body 
Schrodinger equation Hi = — \ T \7 2 + U ext (z). How closely 
the density |^>o(z)| 2 °f the one-body problem approxi- 
mates the density p(z) of the many-body problem de- 
pends on the area density n and the strength of the in- 
teractions, but also on the strength of the trap poten- 
tial. For very tight confinement, the eigenenergies of Hi 
above the ground state energy (or above the first two 
modes in case of a double well) have energies so high 
that their contributions to the iV-body ground state can 
be neglected. The weaker the confinement, the more 
|(/)o(^)| 2 and p(z) will differ from each other. This is 
indeed what we find for the three weaker traps with 
(q, A) = (0.50, 2 • 10 3 ), (0.20, 1 • 10 3 ), and (0.16, 5 • 10 2 ). 
The comparison in Fig. [2] between |0o(z)| 2 and p(z) 
shows that the interactions lead to a wider density p(z), 
that does not agree anymore with the one-body assump- 
tion |</)o(z)| 2 . For the three more confined geometries, 
\4>o i z )\ 2 and p(z) are almost indistinguishable (not shown 
in Fig. [2]). Note that this does not mean that a DDI in 
a tight trap is well described by the one-body Hamilto- 
nian Hi; the opposite is true, in-plane correlations are 
stronger in a tight trap [19]. 

For the layer geometry we define a static structure 
function 5(fc||) as 

) = 1 + J d 3 rd 3 r' e jk " [<?(r, r') - 1] (5) 

= 1 + J dzdz'd 2 r\\ e ik ii r ii \g(z,z',r\\)-l] 

where k\\ is the parallel wave number, g(r,r') is the pair 
distribution function introduced above and kn is any 
wave vector in the xy-planc with wave number k\\ . Note 
that 5(fcy) depends only on fen, while we integrate over 
the z and z' dependence of g(r, r'). In Fig. [3] we show the 
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FIG. 4. (Color online) Integrated two-body density (see text) 
for particles in the same layer pw (bottom panel), and particles 
in different layers P12 {top panel) 



static structure function S(k\\) as a function of fen for the 
six traps studied here. As we go from well-separated thick 
layers to close thin layers we observe a peak in S^fcy) in 
both limits, whereas the peak vanishes in between. It is 
natural to assume that for wide layers the peak is caused 
by correlations due to the intra-layer attraction of the 
dipoles whereas at for a small layer distance it is caused 
by correlations due to the inter-layer attraction of the 
dipoles. However, 5(fey) does not contain enough infor- 
mation to distinguish between the these two mechanisms. 

In order to gain information about intra- and inter- 
layer correlations, we define the partially integrated pair 



5 



densities pn(r\\) and pi2(r\\), 



oo /■oo 



Pn{ r \\) = ^2 / / dzdz' p 2 {r\\,z.z') (6) 

Po Jo Jo 

4 r° r°° 

Pi2{r\\) = ^ / / dzdz' p 2 {r\\,z 1 z') (7) 

Po J~oo JO 

The prefactors are chosen such that pij(fn) — > 1 for 
rii — > oo, thus pn(? , ||) and /O12 (tm ) can be regarded as 
intra- and inter-layer pair distribution functions. They 
are the normalized probabilities to find two dipoles in the 
same layer and in opposite layers at a parallel distance 
ry , respectively, regardless of their z-coordinate within 
the layer. pu(r\\) and pi2{r\\) are shown in the lower and 
upper panel of Fig. [4] for all six traps. For wide, but well- 
separated layers there are strong intra-layer correlations 
at ry =0, whereas the inter-layer correlations are vanish- 
ingly small. This means that two particles in the same 
layer have a very high probability for head-to-tail config- 
urations, with no parallel separation. As we decrease the 
thickness of each layer, these intra-layer correlations van- 
ish. At the same time we decrease the distance between 
layers, thereby increasing the inter-layer correlations. For 
the smallest distance, two particles in different layers are 
strongly correlated and have a high probability for head- 
to-tail configurations. Since the layer is thin, particles in 
the same layer have a vanishing probability for zero par- 
allel separation because of the DDI and the short-range 
repulsion. In both limits of two independent wide traps 
and two close narrow traps the respective strong positive 
correlations of pu and p\2 suggest a tendency towards 
dimer formation, where two dipoles align head-to-tail ei- 
ther within a layer or across two layers. 

What happens, if we would drive the system to even 
larger correlation peaks in pn(r||) or p\2{r\\)7 The insta- 
bility with respect to dimerization manifests itself as a 
numerical instability of the HNC-EL equations. Unlike 
other approximations, the HNC-EL equations have the 
benefit that they do not produce a solution, if a ground 
state of an assumed variational form does not exist. In 
the present case, the Jastrow-Feenberg ansatz Q does 
not allow for the dimerization that our above analysis of 
intra- and inter-layer pair distributions clearly suggests. 
Since the ground state we try to compute does not ex- 
ist, our iterative procedure to solve the HNC-EL equa- 
tions does not converge. In order to actually compute 
the properties of the dimerized phase, one would have 
to optimize a variational ansatz that is flexible enough 
to allow dimerization, or alternatively perform quantum 
Monte Carlo simulations, as e.g. in Ref. [20] . 



III. EXCITATIONS 

A. Bijl-Feynman modes 

Owing to the translational invariance, excitations can 
be characterized by a parallel (i.e. in-plane) wave num- 
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FIG. 5. (Color online) Energy of the first and the second 
excitation mode in Bijl-Feynman approximation for different 
distances of the layers (lines), in comparison to the energies 
obtained for two 2D layers (dots). 



ber fcy . For a given kn there are in principle an infinite 
number of excitations, that are indexed by a perpendic- 
ular quantum number n £ N associated with out-of- 
plane motion. Especially for narrow double-well traps 
these modes have a much higher energy than the low- 
est two modes in the interesting regime of wave numbers 
fcii , therefore we restrict our discussion to the two low- 
est modes and the appearance of a soft mode. In Fig [5] 
we show the first and the second excitation mode, ei(k\\) 
and £2(fc||) in Bijl-Feynman approximation as a function 
of fcii for the six layer geometries for which we studied 
the ground state above. The first and the second ex- 
citation mode are almost degenerate for a large distance 
between layers and considerably split for a small distance. 
Two completely independent DBG layers would of course 
have two-fold degenerate excitation energies. Even for 
the most separated layers, the Feynman dispersion is not 
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FIG. 6. Difference Ae(feii) between the energies of the two 
lowest excitation modes in Bijl-Feynman approximation. 



is completely missing. Positive correlations are possible 
only for the inter-layer pair distribution pn{r\\). The 2D 
results for the two lowest Bijl-Feynman energy disper- 
sions are shown as symbols in Fig. [5j For the wide layers 
that are far apart, the quasi-2D and 2D results differ sub- 
stantially (top left panel), which demonstrates that the 
bending of the dispersion towards forming a roton is an 
intra-layer effect. As we make each layer narrower, the 
quasi-2D and 2D results become almost identical. This 
means that the intra-layer attraction plays less of a role 
and the roton formation is truly an inter-layer effect. 
Note that in the 2D limit, the splitting between the two 
lowest modes vanishes for fen — > in agreement with the 
above argument that the DDI averages out when inte- 
grated over the whole 2D plane. 



truly degenerate, but split for small wave numbers. In 
order to illustrate this, we show the energy difference 
Ae(fcii) = £2(fc||) — ei(A:||) between the two lowest Feyn- 
man energies as a function of k\\ in Fig. [6| for the two 
traps closest to instability. The lifting of the degeneracy 
can only be due to the DDI that is long ranged and hence 
couples even well separated layers. We can estimate the 
typical range of parallel wave numbers for which excita- 
tions are most strongly affected. We assume a circular 
density wave ~ Jo(fcr), of wave number k, in one layer. 
A particle at r = in the other layer will feel a partic- 
ularly strong dipole force if k is such that Jo(fcro) = 
where r is the radius where the DDI changes from attrac- 
tive to repulsive, tq is given by r$ = dtan#, where 9 is 
the angle of the attractive cone of the dipole interaction, 
cos 9 = 1/ V3, which gives r$ = dy2. From this we get an 
estimate for the wave vector k at which we observe the 
strongest dipole coupling, which is k = 2.4048 /(dy/2). If 
we estimate d as the distance between the two maxima of 
the density profiles shown in the top panel of Fig. [T] we 
obtain k s» 4 and k s» 0.24 for the closest and most sep- 
arated layers, respectively. This simple estimate agrees 
reasonably well with the maximum energy splitting of 
the Feynman spectrum at k = 5 and k = 0.35 in Fig. [6] 
Note that for k = the DDI averages out, leading to zero 
DDI-induced splitting for k — > 0, which is what we ob- 
serve for well-separated layers. For the closest layers the 
splitting at k = is large, however, which is caused by our 
short-range repulsion model (er/r) 12 which at such small 
d can be felt between different layers. One could decrease 
a without compromising the stability against intra-layer 
dimerization, but we preferred to tune only the external 
trapping potential while keeping the interaction parame- 
ters fixed. Furthermore, the tunnel splitting is not small 
anymore for the closest layers, adding to the splitting 
caused by the short-range repulsion. 

In order to test our conclusions regarding inter-layer 
coupling for two layers of finite thickness, we also per- 
formed calculations for the limit of two 2D layers. In 
this case the interaction within the same layer is purely 
repulsive ~ r~ 3 and the attractive part of the interaction 



B. Dynamic structure function from CBF-BW 

Calculations of the excitations in the 2D limit of single 
layers have shown [18] that the Feynman approximation 
is adequate for the dispersion relation only at very low 
densities, but correlation effects become more important 
as the density is increased and fluctuations of pair corre- 
lations must be taken into account. Pair correlation fluc- 
tuations are accounted for in the correlated basis func- 
tion - Brillouin-Wigner (CBF-BW) formalism g5J. The 
CBF-BW method not only improves the accuracy of the 
excitation energies, it also describes damping via decay 
of collective modes into lower energy modes. We will 
see that the DDI coupling between layers leads to even 
larger deviations of qualitative features of the excitation 
spectra in the Feynman approximation. 

The CBF-BW method was adapted to layer/film ge- 
ometries in Ref. [SO] and applied to supcrfluid 4 He films 
50 52 and recently to single layers of a DBG .19 . The 
CBF-BW method has been demonstrated to yield excita- 
tion energies much closer to the experimental results than 
the Feynman approximation, even for such a strongly cor- 
related system as 4 He. Further improvement had been 
achieved for bulk 4 He by including fluctuations of triplet 
correlations [S3]. The added complexity, however, pre- 
cludes an application to inhomogeneous systems. 

The CBF-BW excitation energies are conveniently 
obtained by following the linear response approach 
that yields the density-density response operator 
x(r, r',E) and - via the fluctuation-dissipation theo- 
rem |54j - the dynamic structure function S(r, r',i?) = 
— Qtrax(r, r', E)/ir, where E/Ti is the frequency of a small 
external perturbation. The derivation of x{ r , r ',E) can 
be found in Ref. [501 If we project S(r,r',E) onto plane 
waves 

S(k, E) = J d 3 rd 3 r' e lk(r - r,) S(r, r', E) 

we obtain the inelastic cross section for a perturbation 
imparting the momentum Kk to the system. For a given 
k, a peak in S(k, E) at an energy E = E indicates an 
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FIG. 7. (Color online) S(ku , E) is shown for six differ- 
ent traps, with the trap parameters q and A given in each 
panel. For better visibility of low-intensity features, we map 
S(k h E) l/A to a gray scale. The full lines trace undamped 
peaks, the dotted lines are the dispersion relations in Bijl- 
Feynman approximation. 



excitation of energy E. Peaks can have zero linewidth, 
if decay of an excitation is kinematically forbidden, or 
finite linewidth otherwise. Translation invariance in the 
xy-plane implies that the projection ky of k on the xy- 
plane is a good quantum number. A perturbation trans- 
ferring a parallel momentum ftky and energy E to the 
system probes the dispersion relation e n (k\\) of the col- 
lective excitations, that we have calculated above in the 
simpler Feynman approximation. One might think that, 
since only the parallel component of k matters for mea- 
suring the dispersion relation, we can restrict ourselves 
to a parallel k, with a vanishing perpendicular compo- 
nent kj^. However, the corresponding dynamic structure 
function <S(ky , E) probes only excitation modes of even 
symmetry with respect to the rcy-plane. Since we are 
interested not just in the lowest (even) mode but also 
in the second (odd) mode, we will show S(k, E) also for 



wave vectors k which have an angle 9 with the xy-plane. 
A purely perpendicular perturbation (ky = 0) could be 
implemented by fluctuations of the trapping potential ([I]) 
itself, but such a perturbation does not probe the disper- 
sion relation. 



1. Parallel Momentum Transfer- 
In Fig{7] we show S(k, E) for the six different traps 
shown in Fig{TJ the parameters given in table [TI] and 
for wave vectors k that are parallel to the xy-plane, i.e. 
kj_ = and hence fey = k. S(k, E) is represented in Fig. [7j 
by mapping S(k,E) 1 / i to a gray scale. The power of | 
makes sure that also broad, but low-intensity features 
can be seen well. Full lines track peaks of S(k, E) of zero 
linewidth, i.e. which are proportional to a <5-function. 
The resulting line is an undamped dispersion relation. 
For sufficiently large wave number k, the dispersion rela- 
tion merges with the gray area, where damping by decay 
of an excitation into two lower-energy excitations is kine- 
matically possible (i.e. energy and momentum are con- 
served). The Bijl- Feynman spectrum is shown as dotted 
lines for comparison, including also higher modes. For 
wider, well separated traps, the Bijl-Feynman dispersion 
agrees quite well with the CBF-BW result - even for 
the widest trap where a roton starts to form due to the 
intra-layer instability (top left panel). Of course, the 
Bijl-Feynman approximation does not account for damp- 
ing. As we confine the two layers more strongly by in- 
creasing both trap potential parameters A and q, the 
dipole coupling between films leads to a splitting of the 
Bijl-Feynman energies, as discussed above. S(k, E) has 
a much richer structure that is poorly represented by the 
Bijl-Feynman spectrum. On the one hand the density 
increases as the trap tightens (see Fig[T]), and the Bijl- 
Feynman approximation becomes worse at higher den- 
sity. On the other hand, the DDI between layers leads, 
in addition to a splitting of excitation energies, also to 
more decay channels. 

We demonstrate the importance of inter-layer DDI 
coupling by switching it off for the two traps resulting 
in the closest layers (A = 2 • 10 4 and q = 2.0; 2.7). This 
is simply achieved by setting Vdd to zero if z\ and z-i 
have opposite signs. In Figs. [8] and 
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we show £ (fey , E) 

with the full DDI in the left^anels - and without inter- 
layer DDI in the right panels. For q = 2.0 (Fig. [I}, the 
lack of inter-layer DDI almost completely decouples the 
two layers, leading to an almost degenerate Bijl-Feynman 
spectrum. What we get is the dynamic structure func- 
tion of a single layer, which has been studied in Ref. \W\ 
The inter-layer DDI leads to significant additional damp- 
ing for higher energies, seen by the wider peak in the en- 
ergy regime where the dispersion becomes approximately 
quadratic. Note that even without inter-layer DDI, there 
is a bend in the dispersion for q = 2.0, which shows it is 
not so much caused by inter-layer coupling, but by the 
intra-layer repulsion of the DDI that, for much higher 
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FIG. 8. (Color online) The left panels shows S(k\\,E) for 
A = 2 ■ 10 4 and q = 2.0 as in Fig[7] i.e. with inter-layer DDI. 
The right panel shows the corresponding S(k\\ , E) when the 
DDI between the layers is switched off. 



area density, results in the type of roton studied in Ref.lTBl 
in the 2D limit. 




FIG. 9. (Color online) A slice of S(k\\ , E) with (full line) 
and without (dashed line) inter-layer DDI, with A = 2 ■ 10 4 
and q = 2.0. fcii = 5.0 in the upper panel and fcii = 8.0 in the 
lower panel. The vertical lines in the upper panel indicate the 
respective energies of the undamped mode. Arrows show the 
respective excitation energies in Bijl-Feynman approximation. 

While a full S(kt\ , E) map is necessary to track the dis- 
persion relation, the detailed line shapes of the various 
peaks are best seen by plotting slices of S{k\\ , E) for fixed 
values of k\\ ■ The top panel of Fig. [9] shows, for trap pa- 
rameters A = 2 • 10 4 and q = 2.0, a slice of S(k\\,E) at 
fey = 5.0, which is slightly below the value of fey where the 
sharp dispersion curve merges into the damping regime 



and thus becomes broad (see full S(k\\ , E) map in Pig.lsl. 
The full line and dashed line are the results for S(ku 
with and without inter-layer DDI, respectively. The cor- 
responding excitation energies in Bijl-Feynman approx- 
imation are indicated by arrows. The vertical lines are 
the undamped peaks of the sharp dispersion. We see that 
S(k» , E) has only one broadened peak without inter-layer 
DDI, while the inclusion of the inter-layer DDI leads to 
two broadened peaks. We stress again that for parallel 
momentum transfer, S(k,E) only probes the lower, even 
mode, hence the two broad peaks are not due to the split- 
ting of a degenerate eigenmode (5(fcy , E) for non-parallel 
momentum transfer is presented below). 

As fcj| is increased further, the sharp peak loses 
more and more spectral weight and eventually becomes 
damped. This case is shown in the lower panel of Fig. [9j 
where ku = 8.0 and the sharp dispersion has vanished 
for both the coupled and uncoupled bilayer, see Fig. [8] 
Again there is only a single peak without inter-layer DDI 
and two peaks with inter-layer DDI. The lower peak is 
caused by the DDI coupling while the higher one is only 
shifted slightly with respect to its position without DDI 
coupling. Note that the DDI coupling approximately 
doubles the width of the higher peak, hence reduces the 
lifetime of the associated excitation by about a factor of 
two. Thus, as one can expect, the dipole-dipole coupling 
between layers leads to faster decay of excitations com- 
pared to uncoupled layers. 




FIG. 10. (Color online) Same as Fig. |] for A = 2 • 10 4 and 
q = 2.7. 



Finally, in Fig. 10 we compare S(k, E) with and with- 
out inter-layer DDI for even closer layers (A — 2 • 10 4 
and q = 2.7). The bending is now significantly enhanced 
by the inter-layer DDI. In CBF-BW approximation, the 
dispersion (blue line) has a small slope at k\\ = 3, i.e. the 
system is close to "rotonization" . Note that the residual 
splitting of the dispersion without inter-layer DDI is due 
to tunneling and the short-range repulsion as mentioned 
above. 
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2. Non-Parallel Momentum Transfer 

Parallel momentum transfer only probes those exci- 
tations which are even with respect to inversion at the 
z = 0-plane, because a perturbation independent of z is 
even and therefore cannot excite odd modes. In order to 
probe odd modes, we study S^k, E) for wave vectors k 
with an arbitrary angle 9 with respect to the z = 0-plane. 
Fig.[TT]shows 5(k, E) for = 0; 20; 40; 60; 80, for trap pa- 
rameters A = 2 • 10 4 and q = 2.7. We plot S(k, E) as a 
function of fey, not |k|, since only fey is a good quantum 
number that is meaningful for characterizing the excita- 
tion spectrum. Unlike in all previous figures of S(k, E), 
we now add an artificial small imaginary part rj = 0.1 
to the energy E which slightly broadens all features of 
S(k,E). The rationale behind this broadening is that 
it makes the spectral weight of peaks with zero intrinsic 
linewidth visible. 

The case 9 = was shown already in Figj7] (without 
artificial damping), and is shown here again for better 
comparison with 9 > 0. For 6 — indeed only the lowest, 
even mode is visible. As 9 is increased, a second mode be- 
comes visible and gains weight. For 9 = 60, both modes 
can be seen equally clear in 5(k, E) where they appear 
as a narrow dark trace. Note that the second, odd mode 
is damped for small fey. This is very different from the 
low-fen behavior of the lowest mode (sound mode) which 
is not damped because of its negative curvature. The 
damping of the second mode can be seen as a broaden- 
ing (in addition to the artificial broadening) for fen < 1.8. 
We introduce the damping limit E n _ rn {k\\) which is the 
energy above which an excitation of parallel wave num- 
ber fey can decay into two modes with perpendicular wave 
number n and m. E n>m (k\\) is given by 



E„ 



l (k\\) = min[e r 

911 



(fey) +e ro (|q|| 



where, due to the limitations of the CBF-BW approxima- 
tion, e n (Q\\) are the excitation energies in Bijl-Feynman 
approximation, not the excitation energies following from 
CBF-BW itself (inclusion of triplet correlations has 
been shown for homogeneous systems to lead to a self- 
consistent formulation of the self-energy, see Ref. [53]). 
Even modes can decay into combinations where n + m is 
even and vice versa for odd modes. Since we are inter- 
ested only in decays of the lowest two modes, we obtain 
three decay limits, which fulfill £i ; i(fcy) < E 12 (k^) < 
-272,2 (fell)' They are shown in Fig 
lowest mode can decay into (n, m) 
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as dashed lines. The 
(1,1) and (2,2) and 
the second mode can decay into (n,m) — (1, 2). The ef- 
fect of these respective limits are clearly seen in FigfTTj 
Damping indeed sets in as the dispersion relation of the 
mode crosses the damping limit with a symmetry appro- 
priate for the mode. 



Also visible in Fig 11 are interference patterns that 
lead to a modulation of the intensity of S(k, E) as fen 
and thus k± = fey tan# is increased. These are simply 
due to the perpendicular wave number k± being in phase 



or out of phase with even or odd modes. For example, a 
value of fc_L ~ 7r/d, where d is a measure for the distance 
between the layers, leads to a cancellation of the inten- 
sity for even modes, but to a maximal intensity for odd 
modes. 



IV. DISCUSSION AND CONCLUSION 

In this work we generalized our previous studies |19j of 
dipolar Bose gas layers from a single layer in a harmonic 
trap to double-well traps which results in a bilayer ge- 
ometry. As in our previous work, the bilayer calculations 
are based on the HNC-EL method for the many-body 
ground state and on the CBF-BW method for excita- 
tions. Dipolar bilayers have a richer structure than sin- 
gle layers owing to the inter-layer dipole coupling. The 
possibility of head-to-tail pairing of two dipoles on dif- 
ferent layers leads to similar rotonization effects in the 
non-paired (monomer) phase as previously predicted in 
single layers. We restricted ourselves to the calculation 
of ground state correlations of the monomer phase, as 
well as its excitation spectrum, including damping due 
to decay of excitations into two lower excitations. 

We systematically varied the double- well trap parame- 
ters between two close, but thin layers and two well sep- 
arated, but wide layers, while keeping the total area den- 
sity fixed at a modest nr^ = 1. The two end points of the 
range of trap parameters are marked by instabilities of 
the monomer phase. Either if layers are too close or if one 
layer is too wide, inter-layer or intra-layer dimerization 
occurs, respectively. The latter kind of dimers are not 
stable and would quickly collapse via 3-body collisions, 
but inter-layer dimers are stable, given a sufficiently high 
double well barrier. The propensity to pairing was clearly 
seen in the monomer pair distribution functions pi2{r\\) 
or pn (ry), which are the normalized probabilities to find 
two dipoles in different or the same layers, respectively, 
at a parallel distance ry. We showed that, at rii = 0, 
Pi2(f\\) grows a peak for small interlayer distance, while 
pn (ry) grows one if each single layer is sufficiently wide. 
In both cases, the peak of Pij(r\\ = 0) is a precursor to 
the pairing of two dipoles in head-to-tail orientation. 

We presented calculation of the dynamic struc- 
ture function 5(k, E) in the CBF-BW approximation. 
5(k, E) for parallel momentum transfer probes only even 
modes, where we are mostly interested in the lowest one. 
S(k,E) typically consists of a lower, undamped peak 
(that vanishes for higher fey) and two broad peaks that 
are due the inter-layer DDI coupling (without it, there is 
only one broad peak) , which also enhances damping. The 
double-peak structure is not to be confused with the more 
trivial effect that each mode is split into two because the 
inter-layer DDI lifts its degeneracy. The rich structure 
of S{kn , E) is not captured by the simple Bijl-Feynman 
approximation which would predict a single, undamped 
peak for the lowest mode. The intra- and inter-layer in- 
stabilities of the monomer phase are characterized by a 
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FIG. 11. (Color online) S(k, E) is shown as a function of fcii for different angles 9 of k with respect to the plane of the bilayer. 
The trap parameters are A — 2 • 10 4 and q — 2.7 and 5(k, E) was broadened by adding a small imaginary part to the energy, 
r\ = 0.1. As 6 is increased, the second excitation modes becomes visible in 5(k, E). The dots are the energies in Bijl-Feynman 
approximation, and the dashed lines are the damping limits £7 n ,m(fc||) discussed in the text. 



bending of the dispersion relation of the lowest (intra- 
layer dimer) or lowest two (inter-layer dimer) excitation 
modes. This bending, that is less pronounced but still 
visible in the Bijl-Feynman approximation, indicates "ro- 
tonization" , which is well-studied for single layers. As in 
our previous work on single layers, we found that the it- 
erative procedure to solve the nonlinear set of HNC-EL 
equations becomes unstable as we approach rotonization, 
i.e. as the dispersion relation starts to have a local min- 
imum at finite fey . This leads to the conjecture that the 
ground state is only metastable when the excitation spec- 
trum exhibits a roton, while the true ground state, i.e. 
the state of lowest energy is the (intra- or inter-layer) 
dimerized phase. For a proof of this conjecture, however, 
one would need to compare our monomer results with 
results for the dimerized phase to find out which state 
has the lowest energy. Finally, we also presented results 
for non-parallel momentum transfer, i.e. where the angle 
between k and the plane of the layers is non-zero. The 
dynamic structure function depends on both the parallel 
and perpendicular components of k, fcn and k±. k\\ still 
is a good quantum number, while the non-zero k± allows 
to probe also odd modes, particularly the second excita- 
tion mode. Showing S(k, E) as function of the parallel 
wave number, the second mode becomes clearly visible 
for e.g. an angle = 60°. Unlike the lowest mode, the 
second mode is damped for small k» due to decay into 
two excitations of lower energy. 



An interesting topic are the correlations of the 
monomer phase and its excitations generalized to N lay- 
ers. The long-ranged DDI coupling between different lay- 
ers for example lifts an Y-fold degeneracy of the excita- 
tion spectrum and opens many possible decay channels 
for the resulting N modes. Another direction is the study 
of "unbalanced" bilayers where the two layers have differ- 
ent area densities, or bilayers with different kinds of parti- 
cles (e.g. different mass) on each layer. If for example the 
density in one layer is very low, the DDI interaction with 
the other layer would constitute a very well controlled 
model of an impurity particle moving in one layer, cou- 
pled to a bath of particles in the other layer. The mech- 
anism how an impurity attains an effective mass could 
be investigated in a well-controlled fashion over a much 
wider range of densities and interaction strengths than 
in condensed matter systems. 
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